Scenarios of domain pattern formation in a reaction-diffusion system 
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We performed an extensive numerical study of a two-dimensional reaction-diffusion system of 
the activator-inhibitor type in which domain patterns can form. We showed that both multidomain 
and labyrinthine patterns may form spontaneously as a result of Turing instability. In the stable 
homogeneous system with the fast inhibitor one can excite both localized and extended patterns by 
applying a localized stimulus. Depending on the parameters and the excitation level of the system 
stripes, spots, wriggled stripes, or labyrinthine patterns form. The labyrinthine patterns may be 
both connected and disconnected. In the the stable homogeneous system with the slow inhibitor 
one can excite self-replicating spots, breathing patterns, autowaves and turbulence. The parameter 
regions in which different types of patterns are realized are explained on the basis of the asymptotic 
theory of instabilities for patterns with sharp interfaces developed by us in Phys. Rev. E. 53, 3101 
(1996). The dynamics of the patterns observed in our simulations is very similar to that of the 
patterns forming in the ferrocyanide-iodate-sulfite reaction. 
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I. INTRODUCTION 

Pattern formation and self-organization are among the 
most fascinating phenomena in modern science that are 
observed in physical, chemical, and biological systems 
of very different nature (for the books and recent re- 
views on this subject see [p]-p^, where a lot of refer- 
ences to the original works can also be found). As a 
rule, self-organization is associated with Turing instabil- 
ity of the homogeneous state in the nonequilibrium sys- 
tems and spontaneous formation of patterns (dissipative 
structures) in them as the excitation level of the system 
(or some other parameter) is varied [l|-P, p^ , p^ . At the 
same time, when the homogeneous state of such a system 
is stable, by applying a sufficiently strong perturbation 
one can excite static, pulsating, and traveling patterns in- 
cluding solitary patterns — autosolitons (AS) [pyTl|-p^. 

Consider, for example, chemical patterns forming in 
the ferrocyanide-iodate-sulfite (FIS) reaction in a gel re- 
actor [p^5|-pT[ . In a typical experimental setting a pattern 
may form spontaneously as a result of the instability of 
the homogeneous state or can be excited by a short lo- 
calized external perturbation of the system. The pat- 
terns forming in both these cases do not have any qual- 
itative differences. They are essentially the domains of 
high and low concentrations of certain substance sepa- 
rated by relatively sharp walls. The pattern may have 
sophisticated geometry and in general may show very 
complicated spatio-temporal behavior. The properties of 
the patterns do not significantly depend on whether the 
system is monostable or bistable. 

It appears that domain patterns forming in very dif- 



ferent systems may in fact have many common features. 
This has recently been noticed in the case of static do- 
main patterns forming both in equilibrium and nonequi- 
librium systems, where a certain set of domain shapes, 
such as spots, stripes, multidomain and labyrinthine pat- 
terns, and the transitions between them has been ob- 
served ijl^ . The same conclusion can be extended to the 
dynamic patterns in nonequilibrium systems. Indeed, 
traveling, pulsating, self-replicating, and stochastically 
oscillating patterns are observed in the systems as di- 
verse as autocatalytic reactions ||^,|l5|-|l7| , semiconductor 
and gas plasma [pO|-pp^|, or premixed flames [p0|-p2[. 



All this suggests that there exists a universality class of 
the nonequilibrium systems in which pattern formation 
and self-organization scenarios are essentially the same. 

Another important question raised by the experiments 
is to identify the totality of possible types of patterns and 
their behaviors in the systems under consideration, and 
to understand the requirements the system should meet 
in order to be able to produce one type of pattern or the 
other. 

In the present paper we will study a model which is 
a typical representative of the pattern-forming systems 
whose phenomenology was discussed above. We will show 
that the type of patterns that form in this model is deter- 
mined mainly by the relationship between the character- 
istic length and time scales of the systems and the way 
the system is excited. We will also show that by changing 
only these length and time scales and choosing an appro- 
priate form for the external stimulus one can make the 
system form practically all kinds of patterns, both static 
and dynamic, that are observed in the experiments. 
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Our paper is organized as follows: in Sec. II we dis- 
cuss the physical mechanisms of pattern formation phe- 
nomena in reaction-diffusion systems of the activator- 
inhibitor type, using a combustion model as an example, 
and introduce a simple model which we study numeri- 
cally; in Sec. Ill we present the results of a systematic 
numerical study of the reaction-diffusion model and give 
qualitative explanations to the effects seen; in Sec IV we 
use the general asymptotic theory of instabilities devel- 
oped by us in Ref. [g3| to identify the parameter regions 
in which one or the other type of patterns is observed and 
compare these regions with the results of the simulations, 
we also give more substantial quantitative explanations 
for the pattern behaviors that are observed basing on the 
interfacial dynamics approach; and finally, in Sec. V we 
compare our results with various experiments and draw 
conclusions. 



II. THE MODEL 

The model which describes the phenomenology of pat- 
tern formation in many nonequilibrium systems is a pair 
of reaction-diffusion systems of the activator-inhibitor 
type 



Te-g-^=PAd-q{0,r,,A), 



dr) 



(1) 



(2) 



where 9 is the activator, 77 is the inhibitor, I and L are the 
characteristic length scales, and Tg and are the char- 
acteristic time scales of the activator and the inhibitor, 
respectively, q and Q are certain non-linear functions, 
and A is the bifurcation parameter. Equations (|l|) and 
(H) have been extensively used to study pattern forma- 
tion in various nonequilibrium systems. In particular, 
they describe electron-hole and gas plasma, various semi- 
conductor, superconductor, and gas-discharge structures, 
systems with uniformly generated combustion material 
[pO|-p^ ; chemical reactions with autocatalysis and cross- 
catalysis ; models of morphogenesis and population 
dynamics in biology Q . Some systems with phase transi- 
tions, such as diblock copolymer blends and ferroelectric 
semiconductors, are also described by equations which 
can be reduced to Eqs. (0) and (|) 

Pattern formation in the systems under consideration 
is associated with a positive feedback of the activator 
6 which results in "self-production of the activator sub- 
stance", this process of self-production is controlled by 
the inhibitor rj that suppresses the growth of the activa- 
tor. It is these two competing processes that give rise to 
different kinds of patterns in these systems. 

The meaning of the variables 9 and rj can be most eas- 
ily understood for the system with uniformly generated 



combustion material [^T]jT^. Consider combustion pro- 
cess in the flow reactor consisting of a chamber placed 
between the two porous slabs which are being cooled. 
The mixture of fuel and oxidizer is pumped through the 
narrow reactor region between the slabs where it is ig- 
nited. Phcnomenologically, this system may be described 
by the equations for the mass diffusion and the heat con- 
ductance which include the reaction terms, averaged over 
the thickness of the reactor region. If n is the concentra- 
tion of the fuel and T is the temperature of the mixture, 
these equations have the form: 



cp- 



dn 
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DAn + G{n) ~ R{n,T), 



KAT + ER{n, T)-WiT), 



(3) 



(4) 



where G = (no — n)/r„ is the rate of the fuel supply, 
R{n,T) = anexp{—Ea/T) is the reaction rate, Ea is the 
activation energy, E is the reaction heat and a is a coef- 
ficient; W{T) = cpiT — To)/tt is the heat removed from 
the mixture, p is the density, c is the heat capacity of 
the mixture, Tq is the temperature of the environment; 
D and k are the diffusivity and the thermal conductivity, 
respectively, r„ and tt are time constants, and A is the 
two-dimensional Laplacian. In this system the activator 
is the temperature T, and the inhibitor is the concen- 
tration n. Suppose that the system is initially in the 
low-temperature state. Now, if a localized region of the 
mixture is heated up, the rate of the reaction in that 
region will rapidly increase, thus producing more heat 
and igniting the neighboring areas. However, this pro- 
cess cannot go forever, since as the reaction proceeds, 
the fuel is being used up, which in turn decreases the 
rate of the reaction. Thus, a positive feedback is real- 
ized with respect to the temperature and the negative 
feedback with respect to the concentration. If we intro- 
duce the variables 9 = T/Ea, rj = n/no, I = -y/ ktt / cp, 
L — 'JDTn, and A = anoETx/cpEa, we will arrive at 
Eqs. (I) and d). 

Pattern formation is most pronounced in chemical and 
biological systems [0 ||j|Jl^ where the processes of 
self-production of matter are responsible for it. However, 
in the real situation such systems are extremely compli- 
cated. Nevertheless, it is often possible to reduce the de- 
scription of these systems to a pair of reaction-diffusion 
equations of the activator- inhibitor type . A par- 

ticularly simple model of this kind is the "cubic" model, 
which is described by Eqs. (0) and (Q) with 



q = 



V, 



= 9 + 7]- A. 



(5) 
(6) 



This model will be studied numerically and analytically 
in the subsequent sections. 

Kerner and Osipov showed that the properties of the 
patterns and pattern formation scenarios in the systems 
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described by Eqs. (1) and (2) are chiefly determined by 
the parameters e — l/L and a = Tg/rri^ and the shape 
of the nuUchne of Eq. (1) for the activator ||lO|-p^. In 
many cases, including the cubic model and the model 
described by Eqs. (||) and (Q), this nuUclinc is N-shaped 
(see Fig. In these systems static, pulsating, and trav- 
eling patterns may form at different values of the param- 
eters e and a. 

When a ^ 1 and e > 1, that is when the inhibitor is 
slower and shorter-ranged than the activator, only trav- 
eling patterns may exist. In the limit e — *■ oo (or, more 
precisely, for L = 0) the properties of such traveling pat- 
terns were studied in detail in both one-dimensional and 
higher-dimensional cases (see, for example, ||,^,^,|| and 
references therein). In the other limiting case e ^ 1 
and a ^ 1, that is when the inhibitor is long-range and 
fast compared to the activator, only static patterns may 
form ||lO|-p^ . These patterns are essentially the domains 
of high and low activator values separated by the narrow 
walls (interfaces) whose width is of order I. Traveling, 
static, and pulsating domain patterns may form when 
both e <C 1 and a <C 1. 

Let us elucidate the physics of the formation of static, 
traveling, and pulsating patterns, including the simplest 
patterns — AS (Fig. ||), using the combustion system de- 
scribed by Eqs. (|) and (|). First of all, it is clear that if 
both the characteristic time and length scales of the vari- 
ation of the concentration n are much smaller than the 
characteristic scales of the variation of the temperature 
T, no patterns will form. Indeed, if a localized region of 
the system is ignited, all the fuel will burn down very fast 
in that region and the flame will not be able to propa- 
gate, so after a short time after removing the heat source 
the flame will extinguish. Different situation is realized if 
the diffusivity of the fuel is much smaller than the ther- 
mal diffusivity and the characteristic time scale of the 
temperature relaxation tt due to the heat exchange be- 
tween the mixture and the porous slabs is much longer 
than the characteristic time scale t„ of the concentration 
variation determined by the rate of the fuel supply. Then 
the conditions a <C 1 and e 1 may be satisfied and the 
patterns in the form of traveling flames may be excited 
in the system. The existence of traveling patterns is due 
to the fact that because of high heat conductance the 
region of size of order I — ^ ktt jcp around the flame 
is heated up and ignited. As a result, the released heat 
ignites the neighborhood and so on. This leads to the for- 
mation of a flame front moving with the speed v 1/tt 
[Fig. 11(a)] . The fuel after the front burns down, so the 
front is followed by the back separated from the front by 
the distance of order wt„ ^ I. Because of the external 
supply the fuel replenishes at the distances of order vTn 
away from the back. 

If the diffusivity of the fuel is much larger than the 
thermal diffusivity of the mixture, the condition e ^ 1 
and a ^ 1 may be satisfied. Then a traveling flame front 
will stop since high diffusion of the fuel will result in the 
decrease of the fuel concentration ahead of the front and 



lead to the formation of the static pattern in the form 
of a flame cell [Fig. ^(b)]. The existence of such static 
pattern is due to the fact that the flame is maintained 
by the diffusive influx of the fuel from the neighborhood, 
where it is constantly supplied through the porous slabs. 
It is clear that when both e ^ 1 and a ^ 1, pulsating 
flames, or more complex dynamic patterns, may form in 
the system [Fig. §(c)]. The cell may first expand, but as 
it cools down and the fuel is used up in it, the front may 
stop and start traveling back . 

From the physical considerations above follows that in 
order for patterns to be able to form, either e or a has to 
be small. Kerner and Osipov suggested the classification 
scheme for the reaction-diffusion systems of the activator- 
inhibitor type based on the relationship between the val- 
ues of e and a [p|-p^. According to this scheme, a system 
with N-shaped nuUcline of the equation for the activator 
is called f2N system if a <C 1 and e ^ 1; KN system, if 
e <C 1 and a ^ 1; or KfiN systems, if both a <C 1 and 
e ^ 1. Accordingly, only traveling waves (autowaves) 
may form in fJN systems, only static patterns in KN sys- 
tems, and all kinds of patterns in KSIN systems. 

Clearly, the mechanisms of pattern formation in all 
reaction-diffusion systems described by Eqs. (|l|) and (||) 
are essentially the same as those discussed above. For ex- 
ample, similar processes lead to the formation of static, 
traveling, and pulsating patterns in the form of the re- 
gions of high temperature and low concentration of elec- 
trons in the photo-generated electron-hole plasma heated 
in the process of Auger recombination . For this 

reason we may use the simple "cubic" model described 
by Eqs. ^ and (|) with (|) and (H) in all our numerical 
simulations. Depending on a and e, this system will per- 
tain to one of the three classes mentioned above. In this 
paper we will concentrate on the case when the system 
under consideration is either KN or Kf2N system. 



III. SCENARIOS OF PATTERN FORMATION: 
RESULTS OF THE SIMULATIONS 

If we choose L and t,j as the units of length and time, 
respectively, we will write the equations describing the 
cubic model in the form 

= e^Ad -e^' + e + rj, (7) 



^ = A7^-9-7j + A. (8) 

We performed numerical simulations of Eqs. (Q) and 
(^ in two dimensions in a wide range of the parameters e, 
a, and A. As we mentioned earlier, either e or a has to be 
small in order for patterns to be able to form. The case 
of f2N systems was extensively studied by many authors 
l@i^|jlJll ^-iid will not be studied here. We will concentrate 
on the case e <C 1 and arbitrary a. 
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The presence of two very different length scales caused 
by the smallness of the parameter e makes the numeri- 
cal simulations of Eqs. and rather difhcult. The 
simulations were performed on the massively parallel su- 
percomputer. An explicit second order finite difference 
scheme was used to discretize the equations. In order 
to accelerate the algorithm, different grid spacings were 
used for and rj. The boundary conditions were neutral 
or periodic. The grid spacing was chosen in such a way 
that a typical front of a pattern contained about 8 to 
10 points. The decrease of the grid spacing by a factor 
of two resulted in the difference of the distributions of 
9 and 77 by a few percent. No noticeable effects on the 
dynamics were observed. 

Before discussing the results of the simulations, let us 
make a few comments about the system of Eqs. and 
. First of all, it is easy to see that these equations are 
invariant with respect to the transformation 

9^-9, V^-V, A~>~A, (9) 

so one only needs to study the parameter region where 
A < 0. The system under consideration is monostable 
for all values of A. The homogeneous state 9 = 9h and 
77 = rjh with 

0, = 77, = |A|i/3(l - (10) 

is stable for A < Aq — — 1/3\/3 ~ —0.19, where Aq is the 
value of A at which 9fi = Sind rjh — VOi the point on 
the nuUcline of Eq. (^ at which q'g — (see Fig. |l|). It is 
easy to see [p|-[l2[ that for a < 1 the homogeneous state 
becomes unstable with respect to the uniform oscillations 
(Hopf bifurcation) with the frequency 

c^o = ( ) at A>A^ = - ( ) , (11) 

whereas for e < 1 it destabilizes with respect to the fluc- 
tuations with the wave vector fc = fco and zero frequency 
(Turing bifurcation) 

fco = (^7-^) at A>Ac^- • (12) 

Notice that for e — * or a we have Ac ^0 or 
A^ Aq, respectively. Also notice that the homoge- 
neous state is always stable when e > 1 and a > 1. It is 
the considerations of stability of the homogeneous state 
that actually lead Kerner and Osipov to divide the sys- 
tem described by Eqs. (|l|) and (||) into W, KN, and 
KilN systems |l0|-[l2|. Similar classification of pattern- 
forming systems of different nature, including the hydro- 
dynamic systems, was proposed later by Cross and Ho- 
henberg 

Most of our simulations were performed at e = 0.05, 
what can be considered reasonably small. In the first sim- 
ulations we studied the Turing instability. The boundary 



conditions in these simulations were periodic. The initial 
condition was taken in the form of the homogeneous state 
plus small random noise. The system then evolved for a 
considerably long time aX A — —0.1 and a = 0.5 when 
the homogeneous state is unstable with respect to Turing 
instability but stable with respect to the oscillatory insta- 
bility. The process of formation of static Turing pattern 
is shown in Fig. |[ In the early stage {t = 6 and t = 9 in 
Fig. ||) the system nucleates some random distribution of 
the activator and the inhibitor. At the intermediate stage 
{t = 16) the pattern transforms into a number of domains 
with sharp walls, at this point the domains may have ir- 
regular shapes and domain fusion frequently occurs; at 
the late stage {t = 51 and t = 210) the domains rearrange 
so that their shape becomes more regular, and a lot of 
smaller domains die through overcrowding. One can see 
that in the end the static Turing (multidomain) pattern 
consists of many disconnected domains with sharp walls. 
Most of them look like slightly distorted circular domains 
of different sizes, although some are more stripe-like. The 
pattern is metastable, upon longer runs a small portion of 
domains may occasionally disappear or change their ge- 
ometry, but its overall appearance remains the same. The 
interaction between the neighboring domains appears to 
be repulsive, so one should expect an ideal hexagonal 
pattern of circular domains of certain radius to be the 
most stable one. This pattern reminds of the ordered 
cellular flame pattern observed in the combustion exper- 
iments [ pO[ and ordered pattern of current filaments in 
the gas-discharge experiments jl^ . 

In the next simulation the initial condition is taken 
the same, but the parameters now are A = —0.1 and 
a = 0.05, so that the homogeneous state is unstable both 
with respect to the Turing instability and the oscillatory 
instability. The evolution of the system toward the static 
labyrinthine pattern in this case is shown in Fig. ^. The 
early stage of the formation of the pattern {t = 0.6 and 
t = 0.9) is the same as in the previous case, but once 
the domains form they start to oscillate (breathe), so 
the high activator value domains invade almost all the 
space in the time interval from t — 1.2 to t = 2.1, but 
then recede, and the process continues. For this value of 
a the pattern finally stabilizes, and eventually a static 
labyrinthine pattern forms. Notice that both the mul- 
tidomain pattern which formed at t = 955 in Fig. ^ 
and the labyrinthine pattern which formed at t = 9.2 
in Fig. H are perfectly good Turing patterns formed as 
a result of Turing instability, so there is no qualitative 
difference between the two. The form of the Turing pat- 
tern in any particular situation must therefore depend 
on the history of getting the homogeneous system into 
the unstable state. The form of the Turing pattern may 
also strongly depend on the small local inhomogeneities 
always present in real systems p^-p^ . Notice that both 
the multidomain and labyrinthine Turing patterns were 
observed in chemical experiments []6|,p|JT7|. 

Our simulations show that for sufficiently small a the 
uniform self-oscillations of the homogeneous state will al- 
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ways set in upon the destabiHzation of the homogeneous 
state even if the system is unstable with respect to the 
Turing instability. Specifically, for e = 0.05 this will hap- 
pen, if a < 0.02. For larger values of a the static Turing 
pattern will always form, although the oscillations of the 
pattern may last a long time. 

In the next series of the simulations we investigate the 
formation of patterns when the homogeneous state of the 
system is stable. In all these simulations the boundary 
conditions were neutral. 

It is well known that when a <C 1 and \ (fiN sys- 
tems) Eqs. (|l|) and (||) admit solutions only in the form 
of traveling waves (autowaves) PJ^,|5|,|8|, p^ , p^|j2^j29[] , and 
when e 1 and a ^ 1 (KN systems) they admit only 
solutions in the form of static patterns, the simplest of 
which are AS in the form of solitary spots and stripes 
of high or low values of the activator surrounded by the 
"sea" of low or high values of the activator, respectively 
("hot" and "cold" AS) Notice that because of 

the monostability of the considered system the radius of 
a spot or the width of a stripe cannot be greater than cer- 
tain value of order one jl^Jl^. Also note that because 
of the symmetry given by Eq. (|9|), in the system under 
consideration we only need to consider the behavior of 
hot patterns. 

When the homogeneous state of the systems is sta- 
ble, the patterns may be excited by means of sufficiently 
strong external stimulus [ pl]JT^ . According to the gen- 
eral qualitative theory, the excitation level A of the sys- 
tem must be greater than certain threshold value A\, in 
order for AS to be able to form ||lo|-|l^. It is possible to 
show that in the limit e — ^ the value of Ah = —1 in the 
considered system. 

First we will consider KN systems. Our numerical sim- 
ulations show that when A is close to A;,, the initial condi- 
tion in the form of the homogeneous state plus a hot spot 
of size of order several e evolves into an AS in the form of 
a localized static radially-symmetric spot. For e = 0.05 
this happens if Af' < A < A^^ , where A^^ ~ -0.72 

and A).2 ~ —0.55. If the initial condition is taken in the 
form of the homogeneous state plus a hot stripe several e 

^(1) 

where Al"' ~ -0.74 and A)^^ ~ -0.55 (the value of 



wide, it will evolve into a static stripe if A)p < A < A 

~ -0.74 and A^^^ ~ -0.55 (the value of a'^^^ 
obtained from the simulations is rather crude since the 
destabiHzation of the stripe may be incredibly slow). 



If the value of A is increased from to AJ,2^ the ra- 
dius of the spot will grow. However, at certain radius cor- 
responding to the value oi A~ A).2 the spot becomes un- 
stable with respect to the radially-nonsymmetric distor- 
tions of its walls [^12 2^j3^. Qualitatively, this means 
that a "burning spot" tends to ignite the neighboring re- 
gions, but since its radius is bounded from above, the 
only way it can grow in size is by elongation. Precisely 
this phenomenon was observed in the simulations. For 



l(2) 



A < A 



:2) 

c3 



—0.46 the spot elongates and transforms 



become further unstable, leading to the tip splitting and 
the formation of a labyrinthine pattern. This effect is 
seen in the simulation of Fig. H, where A — —0.4, which 
is not much greater than A),^ . There an almost radially- 
symmetric spot ai t — Q elongated and transformed into 
a dumbbell at t = 63 and then further destabilized into a 
more complex shape at < = 98. The process of tip split- 
ting resulted in a more and more complicated pattern at 
t — 162 and t — 292, until the resulting labyrinthine pat- 
tern reached the system boundaries and stopped grow- 
ing Sit t — 985. Similar structures were also observed in 
the simulations in Refs. ||3|,Q and in chemical reactions 
. The labyrinthine pattern that formed in the end 
is a collection of stripe-like domains which are all con- 
nected. Notice that in this simulation a = 0.2. However, 
the smallness of the value of a does not affect the simu- 
lation results as long as a e. On the other hand, this 
choice of a significantly accelerates the simulations. 
In the next simulation we took A = —0.25, which is 

(2) 

further away from A),2 and closer to Ac- The initial con- 
ditions were taken in the form of the homogeneous state 
plus a piece of a curved stripe. One can see (Fig. ^) that 
initially {t = 22) the tips of the pattern grow faster and 
start splitting {t — 87). One can also see that the stripe 
itself start to wriggle, and fingers spring out of the re- 
gions with the highest curvature {t = 118). However, in 
contrast to the previous case, some of the portions of the 
growing labyrinthine pattern detach themselves from the 
body of the pattern. As the time passes, more and more 
portions become detached. In the end the labyrinthine 
pattern that fills the whole system at i = 322 consists 
of 5 disconnected pieces. Notice the great similarity be- 
tween this pattern and the labyrinthine pattern formed 
as a result of Turing instability in Fig. ^. Also notice 
that the disconnected pattern is more likely to form even 

f2) 

for A not very different from A],2 , if a is smaller. This is 
because at small a the dynamics of the pattern becomes 
oscillatory and the labyrinthine pattern may form as a 
result of self-replication of spots, which will be discussed 
below. 

The same picture can be observed, if the initial con- 
ditions are taken in the form of a twisted stripe run- 
ning across the system, if A is big enough. In general, 
for A considerably greater than A[, any localized initial 
condition will lead to the formation of the disconnected 
labyrinthine pattern. However, if these boundary con- 
ditions are used with the values of A closer to A),, a 
wriggled stripe pattern will form in the system. This 
process is related with the fact that a stripe may be- 
come unstable with respect to wriggling of the stripe as 
a whole while being stable with respect to fingering when 



A^2^ <A< A'^^ |1|,|12|,^3U30|. For e = 0.05 the value of 



A^^^ obtained from the simulations is A^^ 



(1) 



-0.32. The 



(2) 

into a stripe. If A > A\^ , the tips of the growing stripe 



evolution of the stripe aX A = —0.45 is shown in Fig. |^. 
One can see that the stripe gets more and more wriggled 
without fingering for a long time. Only when the cur- 
vature of some portion of the stripe becomes sufficiently 
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high, a finger springs out {t = 1418). Notice that finger- 
ing also occurs at the points where the stripe is attached 
to the boundary. At these points the curvature of the 
stripe is high as well. 

Up to now we considered the pattern formation in the 
stable KN system in which the inhibitor is fast. Accord- 
ing to our simulations, the time scale of the inhibitor vari- 
ation does not affect all the results above when a ^ 1 but 
a ^ e. In KfiN systems, in which the inhibitor is slow 
enough, pattern formation scenarios will be qualitatively 
different. Figure || shows the evolution of the system at 
A = —0.4, but a — 0.015 with the initial condition in 
the form of an almost radially-symmetric spot. In con- 
trast to the simulation of Fig. ^ which was performed 
for the same value of A and with the same initial condi- 
tions, instead of transforming into a dumbbell the spot 
splits into two in the course of its evolution {t = 3.3). 
The spots that form split in turn into four {t = 5.2). 
This process of self-replication of spots continues until 
the whole system is filled with the multidomain pattern 
(not shown in the figure), which may stabilize or trans- 
form into a synchronously pulsating (breathing) pattern. 
Notice that self-replication of spots was observed in the 
chemical experiments ^6 17 1 and in the simulations [ pl| . 
A similar phenomenon also seems to occur in the chaotic 
cellular fiames ||2^ . 

Figure |^ shows the value of r] in the center of the sys- 
tem in which the synchronously pulsating multidomain 
pattern formed as a result of spot self-replication. One 
can see that the multidomain pattern forms at relatively 
short times {t ^ 5 for the system 10 x 10), and after that 
the oscillation of the pattern as a whole starts. Some 
restructuring of the pattern occurs later on, what results 
in the changes of the oscillations amplitude and the pat- 
tern's geometry. At times t ^ 120 the pattern's oscilla- 
tions had synchronized and no changes in the oscillations 
amplitude nor in the pattern's geometry were observed 
in the longer runs. 

The mechanism of self-replication can be seen from 
Fig. where a single self-replication event is shown in 
detail. One can see that self-replication is determined by 
the two processes: radially-symmetric pulsations of the 
spot's radius and aperiodic growth of the non-symmetric 
distortion. At the beginning the spot expands as a whole 
(t = 0.6), but at the same time the non-symmetric distor- 
tion builds up {t = 1.4). Then the spot starts to shrink 
in the course of the radially-symmetric pulsations, so the 
connection between its right and left portions gets torn 
a.tt~ 2.1. Att = 2.6 there are two spots looking just like 
the one at t = in the system. Notice that for smaller 
values of A a single self-replication act may take more 
than one pulsation period. 

When the value of a is smaller, the process of self- 
replication of domains may become stochastic, produc- 
ing a kind of turbulence (Fig. p^ ). In the simulation 
of Fig. |l] (a = 0.01,^ = -0.4) the initial condition in 
the form of a small domain initially grows in size, but 
at i = 1.2 local breakdown occurs in its center, so that 



the domain transforms into an annulus. The annulus 
then splits in turn into several smaller domains {t — 5.3) 
which engage into incessant stochastic motion. Each do- 
main is self-replicating, but some of the domains formed 
as a result of this process die as a result of the collisions 
with the other domains, what causes the stochastization. 
Another source of stochastization is the local breakdown 
which occurs, if the domain size becomes too big. The 
number of domains in the system changes randomly with 
time. Each domain is also moving as a whole. The inter- 
action between different domains (and the boundaries) is 
repulsive, so the domain fusion is typically avoided. The 
turbulent pattern is persistent and does not synchronize 
even after long times (Fig. 12). It is observed only at suf- 
ficiently large values of A. When A is relatively small, the 
turbulent pattern usually collapses into the homogeneous 
state after relatively short times. This kind of turbulence 
was observed in the chemical [|l^,0 and combustion |2^] 
experiments. Notice that the turbulence that is observed 
in our simulations is different from the spiral turbulence 
observed by Hagberg and Meron . In our simulations 
we never saw the nucleation of the spiral vortex pairs. 

When the value of a is even smaller, a localized initial 
perturbation transforms into an autowave. In the simu- 
lation of Fig. |l|(a), in which a = 0.007, A = -0.3, the 
domain expands, and at t — 1.1 it transforms into an an- 
nulus which now remains stable and continues to expand. 
This results in an autowave passing through the system 
and annihilating when it reaches the boundaries. In this 
case there is no repulsion between the autowaves. At 
these and smaller a only autowaves form in the system, 
regardless of the value of A. If a random initial condi- 
tion is taken, spiral turbulence typical of the excitable 
autowave media (SIN systems) will form at a ^ e [Fig. 
|l^(b)]. Here the turbulent pattern consists of a random 
arrangement of spiral vortices whose positions are fixed 
in space. The spiral waves always annihilate upon colli- 
sion in this case. This is a well-known phenomenon in the 
17N systems (excitable media), for which L — (3|j|,^. 

Before concluding this Section, let us mention two 
other simulations. In the first [Fig. ^(c)], for which 
a = 0.02 and A = —0.3, a localized initial perturbation 
results in a few replication-like acts, but after that the 
pattern stabilizes into a static disconnected labyrinthine 
pattern. In the second the value of e = 0.2 was taken to 
be not very small (the other parameters are a = 1., and 
A = —0.1), so that the system is away from the asymp- 
totic regime e — > [Fig. |l^(d)]. One can see from Fig. 
|T^(d) that the localized initial perturbation transforms 
into a disconnected labyrinthine pattern in this case as 
well, so, qualitatively, the effects observed in this Section 
are realized when the value of e is not very small. 

The patterns which were described above are the only 
kinds of patterns that were observed in the system under 
consideration. No other types of patterns were observed 
in the simulations, no matter what initial conditions or 
the parameters a, e, and A were used (of course, there are 
"cold" patterns, but in view of Eq. (ph they are equiv- 
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alent to the "hot" patterns studied above). Thus, these 
patterns constitute the totahty of the pattern types of 
the considered system. 



IV. DOMAINS OF EXISTENCE OF DIFFERENT 
TYPES OF PATTERNS AND SCENARIOS OF 
PATTERN FORMATION 

In this Section we will analyze the pattern formation 
scenarios observed in the previous Section, give quan- 
titative explanation for the parameter regions in which 
different patterns form, and explain the transformations 
of one type of pattern to the other on the basis of the gen- 
eral asymptotic theory of instabilities for patterns with 
sharp interfaces developed by us in Ref. , and on the 
basis of the interfacial dynamics approaches developed in 
Refs. [|0|j32[|3|. 

Our simulations of the spontaneous formation of Tur- 
ing patterns confirm the conclusion of the general qual- 
itative theory that at the threshold of Turing instabil- 
ity large-amplitude patterns should form abruptly in the 
system pO| , p^ . According to Eq. (p^), for small e Tur- 
ing instability occurs at A ~ Aq = — 1/3\/3 ~ —0.19, 
with respect to the fluctuations with the wave vector 
kc — eT^I"^. One can see from Fig. || that at early 
stages (t = 16) there are many domains of small size, 
so one could naturally assume that at early stages the 
domain sizes are determined by the wavelength of the 
critical fluctuations, which is of order e^/^. However, as 
can be seen from Fig. |[ at late stages the average size 
of the domains becomes greater and in the end all do- 
mains have roughly the same size. This is not surprising 
since the domains of small size are unstable because of 
the effect of the activator repumping ||l^,|l^ . The thing is 
that because of its long-range character, it is difficult for 
the inhibitor to react on such variations of the activator 
that lead to the expansion of some of the domains and 
the simultaneous shrinkage of their neighbors. This can 
be seen from the estimate of the terms in the dispersion 
relation for the fluctuations around the Turing pattern. 
For simplicity let us consider a hexagonal arrangement of 



circular domains of radius T?.^ with the period Lp. Then 
the fluctuation which leads to the activator repumping 
has the wave vector fc = 7r//!p ^ 1 for £p ^ 1. The 
term in the dispersion relation that causes the instability 
is Ao — —e^/TZl |11 1^,^^, whereas the inhibitor reac- 
tion term, which has the stabilizing effect, is of order e£p 
pO| , p^ , and the values of Cp and TZg are of the same or- 
der. Therefore, when TZg is smaller than TZb ~ e^/^, this 
fluctuation will grow and lead to the expansion of every 
second domain and collapse of the rest, what will result 
in the increase of the pattern's period and the radius of 
the domains. In other words, the domains will grow by 
eating their neighbors until their size and the distance 
between them becomes of order e^/'^ . 

On the other hand, the domain radius cannot be 



greater than TZc2 ~ ^^^^ since at greater radii the domain 
becomes unstable with respect to the non-symmetric de- 
formations and either splits or elongates. The important 
thing, however, is that both TZb and TZc2 are much greater 
than 27r/fcc ~ e^^^ for small e, so the process of forma- 
tion of Turing pattern must always consist of two stages: 
initial domain forming and ripenening. 

Notice that in the presence of small localized inhomo- 
geneities the process of formation of Turing pattern may 
be qualitatively different ||ll|,|l2|. A small localized do- 
main may nucleate at the inhomogeneity, but then as a 
result of the transverse instability of its walls, which oc- 
curs when the domain radius becomes of order e^^^ j23j , it 
will transform into a disconnected labyrinthine pattern, 
if e is not very small, or start to split and replicate itself 
until the system is filled with the domains of size of order 
e^/'^, if e is very small These effects will occur when 
a > e. 

According to Eqs. ( [TTI ) and (|l^), Turing instability 
is the first if a > 2e for small e. However, as we see 
from the simulations, even for smaller values of a, that 
is, when the homogeneous state of the system is unstable 
with respect to both Turing and oscillatory instability, 
static Turing patterns may persist up to smaller values 
of a. For e — 0.05 this happens down to a ~ 0.02. For 
these values of a one can see the competition between the 
Turing patterns and the uniform self-oscillations. For 
a < 0.02 the uniform self-oscillations win, and Turing 
patterns do not form. For a > 0.02 the situation is 
reverse. We did not observe the coexistence of Turing 
patterns and uniform self-oscillations. It is interesting to 
note, however, that the well-formed Turing pattern may 
be stable for even smaller values of a. Also, if there are 
small local inhomogeneities in the system with a ~ e, 
the system will nucleate localized domains before reach- 
ing the instability, as in the case of a ^ e, but then 
the domains will self-replicate and a multidomain pat- 
tern will form in the system. If the value of a is smaller, 
the inhomogeneities will cause nucleation of guiding cen- 
ters |l|,|l|. 

Let us now turn to the patterns that are excited in 
the system with the stable homogeneous state. Our first 
observation is that, as was expected |TT|Jl^ ] and in agree- 
ment with the statements of Fife |Q , any localized initial 
perturbation at first relatively quickly transforms into a 
state in the form of the domain with sharp walls, which 
is the closest to the initial perturbation in shape, and 
then this domain starts to evolve considerably slower ac- 
cording to the equations of the interface dynamics. The 
characteristic time scale for the domain to form is that of 
the activator, that is a, and the characteristic time scale 
of the interface motion is a/e so one can see that as 
long as e ^ 1 the latter time is much longer than the for- 
mer. In this sense one could think that at first the initial 
perturbation evolves into a closest in shape stationary 
state, which then grows into a more complicated stable 
pattern as a result of the instability of that state. In this 
process the early stages of the formation of a complex 
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pattern is determined by the type of the critical fluctua- 
tion with respect to which that stationary state loses its 
stability. Therefore, it is important to know the form of 
possible stationary states and when they become unsta- 
ble. 

The simplest patterns in the considered system are 
static spots and stripes pT| , [T^ , pOt . They are indeed ob- 
served when A is sufficiently close to At,, when a spot- 
like or stripe-like initial perturbations are used, respec- 
tively. We performed numerical simulations of the one- 
dimensional and radially-symmetric versions ofEqs. (|) 
and (^) and found the dependences of the stripe's width 
Cs and the spot's radius TZs versus A at e = 0.05 (Fig. 
|T^ . From these simulations one can see that the solu- 



(1) 



< 



tion in the form of a single static stripe exists at A 
A < Ac, where A^^'' = —0.74, whereas the solution in the 

(2) (2) 

form of a single static spot exists when A, < A < A , , 



where ylf ^ = -0.72 and A"-^"' = -0.24 < A^. When 

(2) 

A^ < A < Ac, the local breakdown occurs in the spot's 
center, so the spot transforms into an annulus. The thing 
is that for e <C 1 the distributions of the activator and 
the inhibitor outside the walls of the spot are related via 
the equation of local coupling |ic|- |l^ , |23|] 



(2) 



q{d,f])^0. 



(13) 



In other words, 9 and rj lie on one of the stable branches 
of the nuUchne of Eq. 6 < 6q in the cold region and 
6 > O'q in the hot region (Fig. |^). As the radius of the 
hot spot grows the value of rj in its center gets smaller, 

(2) 

SO at some value oi A = A)^ < Ac it reaches ?7q, the 
point at which the dependence 9{ri) determined by Eq. 
( p^ ) becomes singular, so a sudden down-jump from one 
branch of the nuUcline to the other occurs, resulting in 
the formation of a new interface in the spot's center and 
the transformation of the spot into an annulus. Notice 
that the process of local breakdown in the center of a spot 
and the formation of an annulus in N systems was stud- 
ied in detail in Refs. |p4] , p^ . Also notice that the same 
mechanism is responsible for the local breakdown in the 
center of a one-dimensional stripe However, it 

does not occur in the particular system we study. 

In higher dimensions spots and stripes undergo insta- 
bilities leading to the growth of certain deformations of 
their walls. Recently we developed a general asymptotic 
theory of the instabilities of domain patterns in arbitrary 
N systems . We have shown that the instabilities are 
determined by the motion and the interaction of the pat- 
tern's walls (interfaces). For sufficiently small e one could 
use the formulas obtained in Ref. |2^] and the depen- 
dences Cs {A) and TZs (A) to determine the critical values 
of A at which one or another instability of the spots and 
stripes occur. The parameters that enter those formulas 
for the considered system are 



B = A, Z = 



2V2 



3 ' ^"2- 



However, for e = 0.05 the agreement between the re- 
sults of the simulations and the predictions of Ref. |2^] 
is rather crude (about 50%). This is due to the fact 
that in the derivations of the critical values of the do- 
main sizes, which are typically of order e^^^, we used 
them as small parameters. However, because of the slow 
e-dependence (1/3 power) this is not a very good assump- 
tion for e ~ 0.05. 

Nevertheless, there is a way to calculate the critical 
values of Cg and TZg which agree with the results of the 
simulations with the accuracy better than 5%. To do 
this, we can use the dispersion relations obtained in Ref. 
p3| in the zeroth order of the perturbation theory in the 
potential V , but not expanding in TZs, or Cs, respectively, 
and keeping the value of C evaluated at A — Ab- The 
latter is because the critical fluctuations are localized in 
the walls of the pattern and this is the way to take into 
account some of the potential V. Having done this, for 
the stripes we have the following dispersion relation 



iauj+e k" + Aq 
eBZ- 



2VCTI2T^ 



[l ± exp{~Cs\/C + k^ + , (15) 



where the constants B, C, and Z are defined in Eq. (|lj), 
k is the wave vector along the stripe, to is the frequency, 
the sign corresponds to the symmetric, the "-" sign 
corresponds to the antisymmetric deformations of the 
stripe, and 



eBjl - e-^-^) 
2ZVC 



(16) 



Similarly, for the spots we will have the following disper- 
sion relation 

9 9 

+ Ao = 



7^? 



-eBZ-^I„,{ns^C + iLu)K„,{nsVC + iuj), (17) 

where Im and Km are the modified Bessel functions, m is 
an integer corresponding to the m-th surface mode, and 
Aq in this case is 

\^ = -^-eBZ-^I,n{TlsVC)Km{TlsVC). (18) 

The instabilities occur when Im w < 0. These tran- 
scendent equations can be solved for Cs and TZs, respec- 
tively, when Im lo — Q, for given k or to. Then, using the 
dependences shown in Fig. |lj, one can find the critical 
values of A at which different types of the instabilities 
occur. 

Let us consider the case of the fast inhibitor a ^ e. 
Then, according to Eq. at A = A^^^ = -0.56 the 

spot will become unstable with respect to the m — 2 
mode, which corresponds to a dumbbell-shaped defor- 
mation. Also, from Eq. (p7| ) follows that the spot desta- 

(14) 

bilizes with respect to the to = mode if A < A, = 
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—0.72. This is in perfect agreement with the results of 
the simulations. 



Similarly, according to Eq. (|15|), the stripe be- 
comes unstable with respect to antisymmetric fluctua- 
tions (wriggling) at A > A^.j'' — —0.61, and with respect 

to symmetric fluctuations (corrugation) at ^ > a''^^ — 
—0.32. The minimum width of the stripe is determined 
by the overlap of the fluctuations of the activator, so it 
is not taken into account in Eq. (|l^) |jll|.|l2|,^ . If we 
do take it into account, we will obtain that the stripe is 
unstable at A < A'"^^'' = —0.74. So, here the agreement 
between the predictions of the theory and the simulations 
is excellent as well. 

The type of the complex pattern that forms in the 
late stages of the destabilization of the simple patterns 
is determined by the dynamics of its interface. For de- 
scription of pattern dynamics in higher-dimensional N 
systems Ohta, Mimura and Kobayashi developed an ap- 
proach which allowed them to reduce the equations sim- 
ilar to Eqs. (U) and (^ to the problem of the interface 
dynamics in the case of slow inhibitor in the limit e — > 
and analyzed the early stages of the transverse instability 
development |Q . Goldstein, Muraki and Petrich derived 
an equation of the interface dynamics for a simple system 
of FitzHugh-Nagumo type in the limit of fast inhibitor 
and weak activator-inhibitor coupling and showed that 
the destabilization of simple patterns lead to the forma- 
tion of the connected labyrinthine patterns |^ . Muratov 
derived the general equation of the interface dynamics for 
N systems described by Eqs. (||) and (||) and showed that 
in the limit e — > and a ^ e only multidomain patters 
must form as a result of the instability and self-replication 
of simple patterns However, because of the slow de- 
pendences of certain parameters on e, multidomain pat- 
terns should in fact form only at e ^ 0.01 in the consid- 
ered system [|3| . Yet the interfacial approach remains a 
good approximation for the dynamics of the pattern for 
e ~ 0.05. In this sense the region 0.01 < e ^ 1 can be 
considered a "crossover" region between the labyrinthine 
and the multidomain patterns. This is the reason why 
we see both connected and disconnected labyrinthine pat- 



la) 
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terns in our simulations. When A is not far from A 
the transverse instability is not very strong, so connected 
labyrinthine patterns form (Fig. Here the stripe 

shape is more favorable than the spot shape. As was 

(2) 

noticed in the previous Section, when A is close to A)^2 
(in the simulations we used A = —0.50), a spot destabi- 
lizes into a single stripe which does not branch. This may 
be qualitatively explained by the following argument. In 
order for branching to occur, a spot has to be unsta- 
ble with respect to the m = 3 mode. According to Eq. 
(O), this should happen at ^ > A^^'' = -0.45. For 

A].2 < A < A).^ the stripe is in turn unstable with re- 
spect to wriggling, so as a result of the instability of a 
spot for those values of A the wriggled stripe (Fig. |^) will 
eventually form in the system. This is precisely what we 



see in our simulations. 

When the value of A gets larger, the transverse insta- 
bility gets stronger and the spot shape becomes more fa- 
vorable. There we see domain splitting predicted in Ref. 
psf and, therefore, the disconnected labyrinthine pattern, 
which is the counterpart of the multidomain pattern in 
this case. We emphasize that this will happen only when 
e is relatively large; according to our simulations, indeed, 
for e < 0.01 only multidomain patterns form in the sys- 
tem. 

Another important thing about the complex domain 
patterns is that the multidomain patterns may coexist 
with the labyrinthine patterns. As was shown by Mura- 
tov, for the same values of the parameters e, a and A 
one could excite both multidomain and labyrinthine pat- 
terns by choosing appropriate initial conditions . For 
example, one could take the pattern that formed in the 
end of the simulation of Fig. |^ and use it as an initial 
conditions for the run with the value of A corresponding 
to the stable homogeneous state. Then in the course of 
the system's evolution the domains will shrink and some 
of them will disappear, but in the end the system will 
be filled with the multidomain pattern similar to the one 
in Fig. H (i = 955), rather than with the labyrinthine 
pattern. 

Recently, Hagberg and Meron studied numerically the 
formation of labyrinthine patterns in a bistable N system 
and explained this effect on the basis of the non-trivial 
properties of solitary fronts that form only in bistable sys- 
tems ||35[|. H owever, according to the experimental obser- 
vations ]17[ and our numerical simulations, labyrinthine 
patterns may form both in monostable and bistable sys- 
tems. In order to make the system of Eqs. (^ and (||) 
bistable one needs to add a coefficient 7 in front of 9 
on the right-hand side of Eq. (^). Then for 7=1 the 
system will be monostable, whereas for smaller 7, for ex- 
ample 7 = 0.5, the system is bistable. We did not see any 
qualitative difference between the patterns in these two 
cases. In the monostable systems the solitary fronts do 
not exist at all, so the domains always have finite width 
at least in one direction. The properties of such pat- 
terns are different from those of the solitary fronts, and 
are essentially determined by the non-local interaction of 
different portions of the pattern's interfaces IIIII2I.I2F 



Besides, in the bistable systems with 7 ^ e and e ^ 1 
the solitary fronts are always unstable with respect to 
the transverse instability. Indeed, according to Eq. ( [T5| ) 
with B — 47, C = 1 + ^7, and Cs = 00, the front 
is unstable with respect to the transverse perturbations 
with the wave vector k ^ e"^/"^. So, in general one can- 
not use the properties of the solitary fronts to explain 
the formation of complex domain patterns in N systems. 
Notice that according to the similar argument, any pat- 
tern whose characteristic size is much greater than e^^^ 
is unstable with respect to transverse perturbations both 
in the monostable and bistable systems. 

Let us now consider the case of slow inhibitor a ^ e. 
Since the prevailing shape in the simulations in this case 
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is a spot, we will look for the instabilities of the circu- 
lar domain. We solved Eq. ([l7| ) in the case a ^ e for 
TO = 0, 1 and 2. We found qualitative agreement with 
the results of the general asymptotic theory for instabil- 
ities of domain patterns ||2^ . 

For TO = the instability occurs at Re uj ^ 0. This in- 
stability leads to the transformation of the static spot 
into a radially-symmetric pulsating (breathing) spot. 
Such pulsating spots were observed in the numerical sim- 
ulations and the experiments [p^|- |l2| , ^6|j3^ ] . The instabil- 
ity occurs when the radius of a spot 7?.™"' < TZg < 71™°"^ , 
where 7?.™*" and T?.™"^ are the functions of a. For 
e = 0.05 and a ^ 0.02 the spot is stable in the whole 
region of its existence. For a < 0.009 the spot is unsta- 
ble for all TZs ■ 

The TO = 1 instability leads to the transformation of a 
static spot into traveling. According to Eq. ([l7|), a spot 
becomes unstable with respect to the to = 1 mode when 
TZs > TZt, where TZt is a function of a and A. Notice 
that the general criterion of such transformations was ob- 
tained by Osipov in Ref. |38|. For the same values of A 
the TO = 1 instability always happens at smaller values of 
a than the to = instability. The instabilities for to > 2 
with Re uj ^ occur at even smaller values of a, when 
the spot is already unstable with respect to to = and 
TO = 1 modes. 

The results of the analysis of the instabilities of simple 
shapes (spots and stripes) and the results of the numer- 
ical simulations can be presented on the diagram (Fig. 
|l5|) . This diagram shows the domains of existence of dif- 
ferent types of patterns in the a — A plane for e — 0.05 
when the homogeneous state of the system is stable. All 
the simulations points, which are marked by Roman let- 
ters in Fig. 15, were obtained by using localized initial 
conditions. The vertical lines in Fig. ^ correspond to 
the values of A at which different instabilities of the do- 
main shapes occur, calculated from Eqs. (|T5| ) and (17). 
One can see that these lines separate the regions in which 
the corresponding types of static patterns are observed in 
the simulations. The letter "s" corresponds to the sim- 
ulations in which the aperiodic relaxation was observed. 
The upper dashed line separates the region in which any 
initial condition relaxes aperiodically to one of the static 
patterns from the region in which the relaxation becomes 
oscillatory. As was expected JlT| , p^ , the transition from 
the aperiodic to the oscillatory relaxation occurs at a ~ e. 
Above the upper dashed line the form of the patterns is 
essentially independent of a; depending on the value of 
A and the initial condition one can see spots, stripes, 
wriggled stripes, multidomain patterns and labyrinthine 
patterns (Figs. |, |, |, 0) . Of course, for e — 0.05 
one should use special initial conditions (not localized) 
to excite the multidomain patterns. 

Below the upper dashed line but above the upper solid 
line the relaxation of the initial excitation of the sys- 
tem results in the formation of a stable static pattern, 
although a few pulsations and spot splittings associated 
with them may occur at the beginning [Fig. |l^(c)], so 



the resulting pattern is disconnected for all values of 

(2) 

A > A),2 . The simulations of this type are marked 
"ps" in Fig. |l^. Notice that in this parameter region 
self- replication of spots does not occur, but the fact that 
the inhibitor is slow makes the domain splitting easier, 
since the inhibitor lags behind the motion of the interface 
driven by the transverse instability of the spot |Q . 

The upper solid line represents the solution of Eq. ( p^ ) 
for TO = and shows the instability line for a spot with 
respect to pulsations (breathing). The simulations show 
that below this line different dynamic patterns form. 
When A is big enough and a is just slightly below the 
upper solid line, spot replication leading to the forma- 
tion of static or synchronously pulsating pattern is ob- 
servedfFig. ||). These simulations are marked by "p" in 
Fig. |l^. As was already mentioned above, for e = 0.05 
self- replication does not occur when the inhibitor is fast. 
Nevertheless, self-replication does occur in the case of 
the slow inhibitor (a ^ e). In this case the transverse in- 
stability, which is the primary cause of the domain split- 
ting and self-replication , is assisted by the instability 
which leads to the radially-symmetric pulsations. This is 
the reason why spot replication does not occur above the 
upper solid line, which is the pulsation instability thresh- 
old for a spot. 

As can be seen from Fig. |l^, the dominant type of dy- 
namic patterns for a ~ e is the turbulence (Fig. |ll|, the 
points marked "t" in Fig. |l^). When the two spots come 
at distances less or of order 1 to each other, the inhibitor 
may not be able to suppress the growth of one spot due 
to the shrinkage of the other (the activator repumping 
effect) what may result in the disappearance of one of 
the spots. However, the surviving spot may self-replicate 
in turn and create another spot. Also, if a spot does 
not have other spots around, it may transform into an 
annulus as a result of the local breakdown in the spot's 
center, and the annulus may then break up into a num- 
ber of spots as a result of the transverse instability. It 
is these three processes uncorrelated in space that make 
the turbulence possible in the KilN systems. 

The turbulence is observed at relatively large values of 
A. This is not surprising. Since the turbulence is caused 
by the self-replication process, it may occur only when 

(2) 

the spot is able to replicate, that is, when A > A).2 . 
For smaller values of A the spots do not self-replicate, 
instead they collapse after a few periods of pulsations 
(simulations marked "c" in Fig. [Tsl ). The meaning of the 
separation between the region where the turbulence and 
the synchronously pulsating patterns are realized is less 
obvious. Qualitatively, the disappearance of some of the 
domains in the course of the pattern's dynamics and the 
local breakdown, the processes that cause the stochasti- 
zation, occur easier when the inhibitor is slower, that is, 
when a is smaller. Of course, in order to make a quan- 
titative explanation of this separation, one has to solve 
a highly nonlinear free-boundary problem in two dimen- 
sions. 
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All the dynamic patterns mentioned above are ob- 
served above the lower solid line, which is the stability 
margin for a spot for m — 1 obtained from Eq. (|l^) . Be- 
low this line a static spot destabilizes and transforms to 
traveling. In this region only autowaves (the simulations 
marked "a" ) form from a localized initial perturbation 
[Fig. |l^(a)]. The autowave patterns that form below the 
lower solid line are essentially the same for all values of a 
(aside from the time and length scales of the pattern) and 
in fact do not differ from the autowaves forming in fiN- 
systems with e ^ 1. This is because at a <C e the diffusive 
precursor is not able to form in front of the traveling pat- 
tern front |jll|,|l2|. Observe that according to Fig. no 
complex static or dynamic patterns (except autowaves) 
form in the system at any A if a < . This fact is in to- 
tal agreement with the general qualitative theory 10 -T^] 
and with the conclusions of the general asymptotic the- 
ory of instabilities . 

Hagberg and Meron explained self-replication of spots 
and formation of turbulence as the consequences of the 
parity-breaking bifurcations [nonequilibrium Ising-Bloch 
(NIB) front transitions] of the planar fronts in bistable N 
systems with the weak activator-inhibitor coupling 
Although their approach is useful for the qualitative or 
heuristic explanation of the formation of the dynamic 
patterns discussed above, it is highly inadequate for mak- 
ing quantitative predictions in general. Indeed, the pro- 
cess of self-replication observed in our simulation cannot 
be viewed as a consequence of local NIB transitions (re- 
versal of the propagation direction of the portions of the 
spot's interface), since, according to our numerical simu- 
lations, the whole spot's interface reverses its propagation 
velocity in the course of self-replication (Fig. |l^). Fur- 
thermore, as can be seen from Fig. ^ the region of the 
system's parameters in which spot self-replication is real- 
ized is determined by the instabilities of the static spot, 
and not of the planar front, or the planar stripe, which 
is the counterpart of the planar front for the monostable 
systems. One can show that according to Eq. ([l5|), both 
the instability of the stripe with respect to pulsations [the 
sign in Eq. (|l5|)] and the instability leading to the 
transformation of the static stripe into traveling [the "-" 
sign in Eq. dla)] lie considerably lower than both solid 
lines in Fig. |15|, which correspond to the respective insta- 
bilities of the spot. Hagberg and Meron predict domain 
splitting and formation of disconnected labyrinthine pat- 
terns only close to the NIB transitions (where a ~ e) 
l^sl , and yet domain splitting and the formation of dis- 
connected labyrinthine pattern occurs solely due to the 
transverse instability far from the presumed NIB tran- 
sitions (Fig. H). Also, as was already mentioned, the 
turbulence that was observed in our simulations is dif- 
ferent from the spiral turbulence observed by Hagberg 
and Meron in the bistable system with relatively weak 
activator-inhibitor coupling [^,^. One could think of 
the turbulence observed by Hagberg and Meron as in- 
termediate between the spiral turbulence observed in the 
excitable autowave media (fiN systems) and the turbu- 



lence observed by us in a KfJN system. In our simula- 
tions the nucleation of spiral vortex pairs is not allowed 
by the local breakdown. The turbulence is produced by 
the constant self-replication of spots with the stochastiza- 
tion caused by the disappearance (annihilation) of some 
of the domains because of their strong interaction with 
the neighbors, and the spontaneous creation of new in- 
terfaces (transformation of a spot into an annulus) due 
to the local breakdown, which occurs, if the size of the 
domain becomes big enough (Fig. |ll|). All this suggests 
that in the general N systems with e <C 1 strong non-local 
interaction between the different portions of the domain 
interfaces and between different domains, high curvature 
of the domain interfaces, the time lag between the mo- 
tion of the interface and the reaction of the inhibitor, 
and the process of local breakdown are crucial and in 
fact determine the type of the pattern that will form in 
the system from a localized stimulus for the given values 
of the system's parameters. 

Before concluding this Section, let us discuss how the 
changes in the values of e should affect the bifurca- 
tion sequences and the pattern formation scenarios in 
the considered systems. As was already mentioned, the 
value e = 0.05 corresponds to a "crossover" between the 
asymptotically small values of e, and the relatively large 
e ~ 1. The value of e = 0.05 is reasonably small to admit 
strong separation of the length scales of the activator and 
the inhibitor, yet it is not very small in the asymptotic 
sense, for which in the considered model we should have 
e < 0.01. It is clear that if the value of e is increased, the 
transverse instabilities will become weaker, so the vertical 
lines corresponding to the instabilities of a spot in Fig. 

will move to the greater values of A. For these values 
of e and a ^ e the stripe shape will become dominant 
over the spot shape, so the typical complex pattern form- 
ing in the system will be the labyrinthine pattern which 
consists of long wrig gled stripes, which may still be dis- 
connected [Fig. 0(d)]. This will also be true for the 
Turing patterns forming as a result of the instability of 
the homogeneous state. If the value of a gets smaller, the 
turbulent patterns will form. Our simulations show that 
in this situation the pattern's oscillations are less likely 
to synchronize than in the case of smaller e because of 
the stronger stochastization due to the local breakdown, 
so it would be easier for the turbulent patterns to form. 
If e is even greater, the patterns can no longer be viewed 
as having sharp interfaces, so the entire phenomenology 
of pattern formation will change. 

On the other hand, if the value of e is decreased, the 
transverse instability will become stronger, and the dom- 
inant shape will become the spot shape. In this case all 
vertical lines which correspond to the instabilities of a 
spot in Fig. 



15 will move toward A = Ah. In the case 



a ^ e only disconnected labyrinthine patterns or mul- 
tidomain patterns will form as a result of the transverse 
instability |33|. The characteristic size of such patterns 
wiU be ei/3^Q^ go at very small e any spot will inter- 
act with many other spots. Because of this interaction, 
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for a e the pulsating patterns will tend to synchro- 
nize more, so the synchronously pulsating multidomain 
patterns should exist in a wider range of the system's 
parameters. Also, as follows from the general asymp- 
totic theory of instabilities |2^, for extremely small e 
(e ^ 10^^) the instability of a spot with respect to the 
fluctuations leading to the formation of a traveling spot 
will always occur at larger values of a than the pulsating 
instability, so it is natural to conclude that complex dy- 
namic patterns, such as pulsating multidomain patterns 
and turbulent patterns, cannot be excited by a localized 
stimulus at such small values of e. 



V. CONCLUSION 

In this paper we performed a complete numerical study 
of different types of domain patterns in a two-dimensional 
N system and investigated all major scenarios of their 
formation. We confirmed the conclusions of Kerner and 
Osipov [lC[-|l^ that the type of the domain patterns 
forming in the N systems is determined mainly by the 
two basic parameters of the system: a = T^jr^ and 
e = l/L = DgTo / DrfTri- Thcsc parameters are the ra- 
tios of the characteristic time and length scales of the 
variation of the activator and the inhibitor and are de- 
termined by the local kinetic coefficients: the relaxation 
times and the diffusion coefficients Dg and Drf. In real 
systems these kinetic coefficients strongly depend on the 
excitation level of the system and the state of the en- 
vironment, the presence of small amounts of impurities 
or catalysts which, for example, may change the rates of 
recombination of nonequilibrium carriers in semiconduc- 
tors or the rates of chemical reactions, and so on. For 
example, by varying only the temperature of the semi- 
conductor lattice, one can significantly change the criti- 
cal parameters of the electron-hole plasma described by 
Eqs. d) and d) H. 

We have shown that by changing the values of e and 
a and the control parameter A (in the physical systems, 
such as electron-hole plasma, A is the system's excitation 
level), that is, in essentially the same system the whole 
variety of domain patterns and pattern formation scenar- 
ios is realized. This general conclusion explains theoret- 
ically the results of recent experiments by Lee et al. on 
the FIS reaction [p^^ where they showed that in this 
chemical system by changing relatively weakly its chem- 
ical composition and the form of the initial perturbation 
it is possible to excite all major types of domain patterns 
and see various scenarios of their formation, which are 
qualitatively the same as those observed in our simula- 
tions. 

In KN systems, that is, when e ^ 1 and a ^ e only 
static patterns form. When the homogeneous state of 
the KN system is stable, these patterns can be excited 
by applying a sufficiently strong perturbation (hard exci- 
tation). We found that by changing only the parameter 



A and the form of the initial perturbation one can excite 
localized spots and stripes, connected (Fig. ^ and dis- 
connected (Fig. |6|) labyrinthine patterns, wriggled stripes 
(Fig. ^ and multidomain patterns (Fig. All these 
patterns were found in the chemical experiments . 
Ordered multidomain patterns were also observed in the 
high frequency gas-discharge experiments | [l9| ] and in the 
combustion experiments [Q. We also found that at the 
same parameters of the KN system it is possible to ex- 
cite a great variety of shapes of the static patterns by 
changing only the form of the initial perturbation. In 
particular, it is possible to excite the labyrinthine and 
multidomain patterns in the system with the same val- 
ues of the parameters e, a, and A ||3^. This variety of 
domain shapes is observed when e is not very small, since 
for very small values of e only disconnected multidomain 
patterns will form [^ . However, such small values of e 
can hardly be realized in a typical experimental situation. 
That's why in our paper we paid particular attention to 
the case e — 0.05. 

Static domain patterns in KN systems may also form 
spontaneously as a result of Turing instability of the ho- 
mogeneous state. In the ideally homogeneous KN sys- 
tems these patterns are as a rule quasiperiodic (Fig. ||) 
and in general their period is not the same as the period 
of the critical fluctuation 27r/A:o [see Eq. (12)], but is de- 
termined by the stability of the pattern. In real systems 
the type of the pattern will be determined by the small 
local inhomogeneities. Depending on the form of the in- 
homogeneity, and also on the system's parameters, all 
types of static domain patterns will form spontaneously 
in KN systems. 

In riN systems, that is when a <C 1 and a < e^, 
only the uniform relaxation self-oscillations may form 
spontaneously. In such systems with the stable homo- 
geneous state various autowave patterns, including ex- 
panding traveling waves [Fig. [l3|(a)] and spiral waves, 
can be excited by an external perturbation. If the initial 
conditions are sufficiently random, the spiral turbulence 
[Fig. ^(b)] forms in VtN systems. These patterns were 
first discovered by Zaikin and Zhabotinsky in an oscilla- 
tory chemical reaction [^^ and have been subsequently 
studied for more than two decades in a variety of systems 
[^,||j|, || . Notice that the autowave patterns are also ob- 
served in the FIS reaction [|l^ . 

The most diverse picture of pattern formation is ob- 
served in the case of KfiN systems, that is, when e ^ 1 
and < a < e. In these systems both the uniform 
self-oscillations and the Turing patterns may form spon- 
taneously. In the KfiN systems with the stable homoge- 
neous state it is possible to excite static, traveling, and 
pulsating (breathing) patterns. We showed that in these 
systems the remarkable effect of self-replication of spots 
(Figs. H an d [l0| ) recently discovered in the same FIS 
reaction [ p6| , [l7[ is realized. Depending on the system's 
parameters, this process may lead to the formation of 
static or pulsating multidomain pattern, or to the forma- 
tion of turbulence which is qualitatively different from 
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the spiral turbulence observed in systems. This tur- 
bulence consists of random creation and annihilation of 
spots. Precisely this kind of turbulence was observed in 
the experiments in the very same FIS reaction ||l7|] , and 
also in the combustion experiments . 

We would like to emphasize that the whole variety of 
pattern formation scenarios observed in our numerical 
simulations and illustrated in Fig. |l^ is explained with 
the remarkable accuracy by the asymptotic theory of in- 
stabilities of domain patterns in reaction-diffusion sys- 
tems This is because the scenarios of pattern for- 
mation are determined by the two processes: at first any 
initial perturbation quickly transforms to a state close to 
some stationary state, then the evolution of that quasis- 
tationary state is determined by its instabilities and the 
form of the critical fluctuations with respect to which 
the instability is realized. Indeed, the stationary states 
correspond to the solutions of Eqs. and (||) with the 
right-hand sides equal to zero, so when we use the initial 
conditions which are significantly different from any sta- 
tionary state, we in fact make the right-hand sides of Eqs. 
(|I|) and large and, therefore, the time derivatives of 6 
and J]. The large time derivatives will cause such changes 
in 9 and 77 which will lead to the transformation of the ini- 
tial condition to the state close to some stationary state, 
in which the time derivatives of and rj are small. For 
example, if a square of relatively small size is used as the 
initial condition, it will first transform to a state close to 
a spot; if a square of large size (3> L) is used, it will first 
transform into an annulus as a result of the local break- 
down in the center of the square (this effect was studied 
in detail in Refs. |]l0|-p^). If the system is bistable, the 
initial condition in the form of a large square may trig- 
ger the wave of switching from one stable homogeneous 
state to the other. The evolution of these states close the 
stationary states, and, therefore, the pattern formation 
scenarios, will be determined by their stability: a spot 
may form from a square of small size if the parameters of 
the system are such that it is stable, or an annulus may 
form from a square of large size. The radius of the spot, 
or the width and size of the annulus will be determined 
by the corresponding stable states and will only depend 
on the parameters of the system (such as e and A). If 
the parameters of the system are such that these states 
happen to be unstable, then, depending on the type of 
the instability, which is determined by the parameters e, 
a, and A, all kinds of complex patterns will form. 

One of the bright pattern formation scenarios consists 
of the transformation of the localized excitations into the 
patterns that fill the entire system. In certain sense one 
could think of this as the self-completion of a pattern 
from its small fragment. According to Figs. H] - the 
form of such patterns depends on both the initial condi- 
tions and the integral parameters of the system (in our 
case e, a, and A) which are determined by the system's 
kinetics. It is important that by changing only these 
integral parameters it is possible to excite qualitatively 
different patterns for the same initial condition. Thus, 



the systems we consider have a remarkable prope rty — 
a kind of associative memory (see also Refs. |ll| , [l2| ): the 
form of a pattern is determined by the integral param- 
eters of the ideally homogeneous system, and they can 
be reconstructed with certain probability from a small 
fragment (sufficiently localized initial perturbation). 

Patterns of the same morphology as those in the KN 
system studied by us also form in a variety of the equi- 
librium systems, such as garnet ferromagnets, ferroelec- 
tric and ferrofluid films, Langmuir monolayers, phase- 
separating copolymer blends (for a recent review and the 
references see Ref. @|). Our simulations suggest that 
the complex domain patterns, such as multidomain or 
labyrinthine patterns, are driven by the dynamics of their 
interfaces coupled to the long-ranged inhibitor field. It is, 
therefore, natural to expect qualitatively the same pat- 
tern formation scenarios in the equilibrium systems with 
the competing repulsive interactions and the strong sep- 
aration of length scales. In the nonequilibrium systems 
the inhibitor does not necessarily react on the motion 
of the pattern interfaces instantaneously. The time lag 
of the inhibitor makes the existence of complex dynamic 
patterns possible in the nonequilibrium systems. 

The only kinds of patterns that form in the exper- 
iments with the cellular fiames |^l| and with the gas- 
discharge system we did not see in our numerical 
simulations are traveling spots and hopping patterns. Re- 
cent work of Krischer and Mikhailov suggested that a 
sufficiently strong global coupling, which is absent in our 
model, might be needed to see these kinds of patterns 
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FIG. 2. The simplest domain patterns: (a) traveling, (b) 
static, (c) pulsating one-dimensional autosolitons. 

FIG. 3. Formation of static Turing pattern. The dis- 
tributions of the activator at times. The parameters used: 
£ = 0.05, a = 0.5, A = -0.1. The system's size is 20 x 20. 

FIG. 4. Formation of static Turing pattern in an oscilla- 
tory system. The distributions of the activator at different 
times. The parameters used: e = 0.05, a = 0.05, A = —0.1. 
The system's size is 20 x 20. 



FIG. 5. Formation of a connected labyrinthine pattern. 
The distributions of the activator at different times. The pa- 
rameters used: e = 0.05, a = 0.2, A = —0.4. The system's 
size is 20 x 20. 



FIG. 6. Formation of a disconnected labyrinthine pattern. 
The distributions of the activator at different times. The pa- 
rameters used: e = 0.05, a = 0.2, A = —0.25. The system's 
size is 20 X 20. 



FIG. 7. Formation of a wriggled stripe. The distribu- 
tions of the activator at different times. The parameters used: 
e = 0.05, a = 0.2, A = -0.45. The system's size is 20 x 20. 

FIG. 8. Self-replicating spots. The distributions of 
the activator at different times. The parameters used: 
e = 0.05, a = 0.015, A = -0.4. The system's size is 20 x 20. 

FIG. 9. The value of rj in the center of the 

system as a function of time in the simulation with 
e — 0.05, a = 0.015, A — —0.3 which resulted in the pulsating 
multidomain pattern. The system's size is 10 x 10. 

FIG. 10. A close-up of a self-replicating spot. Same sim- 
ulation as in Fig. (^). The region shown is 8 x 8. 

FIG. 11. The onset of turbulence. The distributions 
of the activator at different times. The parameters used: 
e = 0.05, a = 0.01, A = -0.4. The system's size is 10 x 10. 

FIG. 12. The value of rj in the center of the system as a 
function of time for the simulation of Fig. [lil. 



FIG. 1. The nuUclines of Eqs. (||) and (|). 
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FIG. 13. The distributions of the activator at different 
times for different processes: (a) formation of an autowave 
(e = 0.05, a = 0.007, A = -0.3, the system is 20 x 20); (b) for- 
mation of spiral turbulence (e = 0.25, a — 0.02, A = —0.3, the 
system is 100 x 100); (c) stabilization of a pattern formed in 
the process of splitting of spots (e — 0.05, a = 0.02, A = —0.3 
the system is 10 x 10); (d) formation of a complex pattern 
outside the asymptotic regime (e = 0.2, a = 1.0,^ = —0.15, 
the system's size is 40 x 40). 



FIG. 14. The dependences of Ca on A for the 

one-dimensional AS (a) and TZs on A for the radi- 
ally-symmetric AS (b) for e — 0.05. Results of the numerical 
simulations of Eqs. (0) and (^). 



FIG. 15. The domains of existence of different patterns 
at e = 0.05. The upper part of the figure shows the regions 
where the corresponding patterns may exist at a S> e. The 
lower-case letters indicate the long-time behavior of the sys- 
tem at the parameters corresponding to the position of the 
letter when the system is locally excited at t — 0: "s" — 
the system aperiodically relaxes to a static pattern; "ps" — 
the system relaxes to a static pattern, but the relaxation has 
oscillating character, a few splittings may occur at the begin- 
ning; "p" — as a result of self-replications at the beginning 
a stationary pulsating (breathing) pattern forms in the sys- 
tem; "c" — the initial excitation collapses as a result of the 
growing amplitude of pulsations; "t" — turbulence develops 
in the system; "a" — the initial excitation transforms into 
an autowave traveling through the system and disappearing 
at the boundaries. The upper dashed line shows schemati- 
cally the region where the character of pattern's relaxation 
changes from aperiodic to oscillating. The upper solid line is 
the threshold of the instability with respect to the uniform 
pulsations for a radially symmetric AS, obtained from Eq. 
(p^. The lower solid line is the threshold of the instability 
leading to the transformation of the radially symmetric AS 
into traveling, obtained from Eq. (|l7|). The lower dashed 
lines show schematically the borders of the parameter regions 
where pulsating patterns or turbulence are realized. 
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